%----------------------------------------------------------------------------------------
%	PROBLEM 1
%----------------------------------------------------------------------------------------

% To have just one problem per page, simply put a \clearpage after each problem
\newpage
\begin{homeworkProblem}
\section{Simulated Annealing} \label{sa}
\subsection{Problem statement}
Implement two stochastic local search (SLS) algorithms for the traveling salesman problem with time windows (TPSTW), building on top of the perturbative local search methods from the first implementation exercise.
\begin{enumerate}
  \item Run each algorithm 25 times with different random seed on each instance. Instances will be available from http://iridia.ulb.ac.be/˜stuetzle/Teaching/HO/. As termination criterion, for each instance, use the maximum computation time it takes to run a full VND (implemented in the previous exercise) on the same instance and then multiply this time by 1000 (to allow for long enough runs of the SLS algorithms).
 \item Compute the following statistics for each of the two SLS algorithms and each instance:
 \begin{itemize}
   \item Percentage of runs with constraint violations
   \item Mean penalized relative percentage deviation
 \end{itemize}

\item Produce box-plots of penalized relative percentage deviation.
\item Determine, using statistical tests (in this case, the Wilcoxon test), whether there is a statistically significant difference between the quality of the solutions generated by the two algorithms.
\item Measure, for each of the implemented algorithms on 5 instances, the run-time distributions to reach sufficiently high quality solutions (e.g. best-known solutions available at http://iridia.ulb.ac.be/˜manuel/tsptw-instances\#instances).
Measure the run-time distributions across 25 repetitions using a cut-off time of 10 times the termination criterion above.
\item Produce a written report on the implementation exercise:
\begin{itemize}
  \item Please make sure that each implemented SLS algorithm is appropriately described and that the computational results are carefully interpreted. Justify also the choice of the parameter settings and the choice
of the iterative improvement algorithm for the hybrid SLS algorithm.
  \item Present the results as in the previous implementation exercise (tables, box-plots, statistical tests).
  \item Present graphically the results of the analysis of the run-time distributions.
  \item Interpret appropriately the results and make conclusions on the relative performance of the algorithms across all the benchmark instances studied.
\end{itemize}
\end{enumerate}


\subsection{Introduction} \label{sec:introSA}
The simulated annealing is a simple SLS method named after the homonym physical process used in material science to obtain desired properties on the processed material.

This process involves heating material to above its critical temperature, then supervise its cooling by means of a precise schedule in order to control the temperature and thermodynamic free energy.

The variations in those two quantities across the different states of the physical system will determine the structure and the properties of the material at the end of the cooling.

The simulated annealing method consider the different candidates solutions as the state of the system, and tries to reach a ground state (i.e. a globally optimum solution) through transitions that are controlled by the virtual temperature ($T$ parameter) and the virtual thermodynamic free energy of the current state ($f(s)$).

The general outline of the implemented algorithm is the following: 

\begin{algorithm}[!h]
  \caption{Simulated Annealing - Outline}\label{sa:outline}
  \begin{algorithmic}[1]
    \State \emph{GenerateInitialSolution}
    \State \emph{InitalizeTemperature} 
    \While{!(\emph{TerminationCondition})}
        \State $s' \gets $\emph{ProposalMechanism}($s$)
        \If{\emph{AcceptanceCriterion}($s,s',T$)}
          \State $s \gets s'$
        \EndIf
        \State \emph{UpdateTemperature}
    \EndWhile
\end{algorithmic}
\end{algorithm}

As in the physical process, the most important parameter in the process is the cooling schedule, which in this case will  directly determine the temperature variation across iterations, and indirectly the quality of the current solution.

In fact the temperature parameter is used in the probabilistic acceptance criterion which will determine the state transition, and will influence the probability of accepting worsening solution throughout the process.  

\subsection{Algorithm structure} \label{sec:algstrucSA}
The proposed algorithm is an implementation of the general Simulated Annealing metaheuristic (cf. \cite{kirkpatrick1983optimization}), taking some design choices especially tailored for the TSPTW problem from \cite{ohlmann2007compressed}.
The main design choices of the algorithm are the following:
\begin{itemize}
  \item \emph{GenerateInitialSolution}: Random generation of a permutation of nodes.
  \item \emph{InitalizeTemperature}: Parameter of the simulation.
  \item \emph{TerminationCondition}: As in \cite{ohlmann2007compressed} with time bounds.
  \item \emph{ProposalMechanism}: Random uniform choice in the 1-shift insert neighborhood.
  \item \emph{AcceptanceCriterion}: Metropolis condition considering as improving also a shorter but infeasible tour.
  \item \emph{TemperatureUpdate}: As in \cite{ohlmann2007compressed}.
\end{itemize}

The design choices can be justified as follows:

\begin{itemize}
  \item \emph{GenerateInitialSolution}: The heuristic developed in the first implementation exercise construct a solution by applying random perturbations to the tour generated by ordering the nodes according to their time windows' closing time.
  
  In other words, the starting point of the algorithm is always in the neighborhood of the solution defined by the ordered permutation.
  
  This introduces a strong bias in the solution research, which may lead the algorithm to converge to sub-optimal solutions.
  
  To remove this bias, the generation of an initial solution is made by generating a random permutation of nodes. 
%   \item \emph{InitalizeTemperature}: The initial value of the temperature is determined by analyzing the local topology of the search space, using techinques described in the literature (cf. \cite{ohlmann2007compressed}).
   \item \emph{TerminationCondition}: The termination condition sets a lower bound on the number of computations in order to obtain a solution of a sufficiently high quality.
  Since the best solution is known and the run-time distribution analysis is required, the algorithm execution is stopped once the maxmium run-time is reached.
  \item \emph{ProposalMechanism}: Empirical tests on the given instances has shown that the 1-shift neighborhood allows to reach higher quality solution with respect to the 2-swap (transpose or exchange) one.
  \item \emph{AcceptanceCriterion}: The acceptance criterion is a modified version of the original metropolis condition in order to take into account the time windows constraints.
  \item \emph{TemperatureUpdate}: Geometric temperature decrease.
\end{itemize}


\begin{algorithm}
\caption{Simulated Annealing TSPTW}
\label{SA:TSPTW}
\begin{algorithmic}[!h]
\Procedure{SimulatedAnnealing}{$\alpha,T_{0},l,t_{\max},f_{best},$}    
  \Require Annealing schedule
  \State $s \gets$ \emph{RandomPermutation}()
  \State $T \gets T_0$  
  \While{!\emph{TerminationCondition}($s,t_{\max},f_{best}$)}
    \State $s' \gets$ \emph{ProposalMechanism()} 
    \If{\emph{UniformRandom}() $<$ \emph{AcceptanceCriterion}($s,s',T$)} 
      \State $s \gets s'$
    \EndIf
  \State $T \gets$ \emph{Update}($T$) 
\EndWhile
  \State \textbf{return} $s$
\EndProcedure
\end{algorithmic}
\end{algorithm}

\begin{center}
  
\begin{minipage}{.45\textwidth}
\centering
\begin{tikzpicture}[->,>=stealth',shorten >=1pt,auto,node distance=3cm,
                    semithick]
  \tikzstyle{every state}=[fill=blue!20,draw=none,thick]

  \node[initial,state]  (IS)                    {\emph{IS}};
  \node[state]          (SA) [right of=IS]      {\emph{SA}};
  \node[state,accepting](End)[below left of=SA] {$\Omega$};

  \path (IS) edge              node {DET:IT} (SA)
        (SA) edge [loop above] node {CDET(not \emph{TC})} (SA)
             edge              node {CDET(\emph{TC})} (End);
\end{tikzpicture}
\end{minipage}%
\hspace{1.5cm}
\begin{minipage}{.45\textwidth}
\centering
\paragraph{Nodes}
\begin{itemize}
  \item \emph{IS} $\equiv$ \emph{InitialSolution}()
  \item \emph{SA} $\equiv$ \emph{SimulatedAnnealingLoop}
\end{itemize}
\paragraph{Conditions}
\begin{itemize}
  \item \emph{TC} $\equiv$ \emph{TerminationCondition}($s,t_{max},f_{best}$)
  \item \emph{IT} $\equiv$ \emph{InitializeTemperature}()
\end{itemize}
\end{minipage}
\captionof{figure}{Simulated Annealing GLSM}
\end{center}


% \subsubsection{Initial solution}
% \begin{algorithm}
% \caption{Initial solution}
% \label{SA:Init}
% \begin{algorithmic}
% \Procedure{Initial Solution}{}

% \State $p \gets$ \emph{UniformRandom}($1,\frac{|N|}{10}$) \Comment{$p$ : Number of perturbations}
% \State $s \gets$ \emph{SortCitiesUBTW}() \Comment{Sort nodes based on upper bound of time window}	

% \While{$p > 0$}
%   \State $p_p \gets$ \emph{UniformRandom}($1,|N|-1$) \Comment{$p_p$ : Perturbation point}
%   \State $p_i \gets$ \emph{UniformRandom}($1,p$) \Comment{$p_i$ : Perturbation intensity}
%   \State $s \gets$ \emph{Shuffle}($1+p_p,1+p_p+p_i$) \Comment{Shuffle the solution components occuping the position bounded by the indexes given as parameters}
%   \State $p \gets p-1$
% \EndWhile	

% \State \textbf{return} $s$

% \EndProcedure    
% \end{algorithmic}
% \end{algorithm}

% The initial solution is generated using the heuristic developed for the previous implementation exercise.
% The tour is first constructed by ordering the cities according to their time windows' closing time, then a stochastic perturbation is applied.
% The modification consist of randomly selecting a number $p$ of perturbation points in the tour, then shuffling the order of the neighboring cities of those point, each time considering a randomly chosen number of neighbors.
% In this way, a solution is built by taking into account the heuristic information (to minimize the number of constraint violations) as well as a stochastic component that allows to differentiate the initial starting point across different runs of the algorithm.

% \subsubsection{Temperature initalization}
% \begin{algorithm}
% \caption{Temperature initialization}\label{saTInit}
% \begin{algorithmic}
% \Procedure{InitializeTemperature}{$\chi_0$}
%   \State $a \gets 0$
%   \While($i < r$)
%     \State $s_1 \gets$ \emph{RandomPermutation}()
%     \State $s_2 \gets$ \emph{ProposalMechanism}()
%     \State $a \gets |f(s_1)-f(s_2)|$
%     \State $i \gets i + 1$
%   \EndWhile
%   \State $\bar{f(s)} \gets \frac{a}{i+1}$
%   \State \textbf{return} $\frac{\bar{f(s)}}{\ln(\frac{1}{\chi_0})}$
% \EndProcedure
% \end{algorithmic}
% \end{algorithm}

% The initial temperature value is determined by taking into account the local topology of the current instance, as well as the initial rate of uphill transition (i.e worsening moves) accepted ($\chi_0$).
% An estimate of the local topology is made by drawing randomly $r=5000$ solutions and computing the average difference among the evaluation of these solutions and their 1-shift neighbors.
% The initial temperature value is then computed according to the formula described in \cite{ohlmann2007compressed}, as the return instruction shows.

\subsubsection{Termination condition}
\begin{algorithm}[H]
  \caption{Termination Condition}\label{SA:Term}
  \begin{algorithmic}[1]
    \Procedure{TerminationCondition}{$s,t_{\max},s_{best}$}
       \If{ $(f(s) = f_{best} \wedge \Omega(s) = 0)$}
				      \State \textbf{return true}
			 \EndIf
			  
			  \If{ $t > t_{\max} $}
				      \State \textbf{return true}
			  \EndIf
       
       \If{$T_{u} > T_{u_{\min}}$}
			      \If{ $T_{us} > c_{t} $}
				      \State \textbf{return true}
			      \EndIf
			\Else    
        \State \textbf{return false}
      \EndIf
    \EndProcedure
\end{algorithmic}
\end{algorithm}

The termination condition of the algorithm is an adaptation of the one proposed in \cite{ohlmann2007compressed}, including also the discovery of the global optimum and the exceeding of the maximum computation time $t_{\max}$ as stopping condition.
The original termination condition consists of number of temperature updates $c_t$ that must be performed without observing any changes in the current best tour.
It also included a lower bound on the algorithm run-time by fixing a minimum number of temperature updates $T_{u_{\min}}$ that have to be performed before stopping the algorithm.
In this implementation, the maximum computation time constraint has an higher priority with respect to the original termination condition.

\subsubsection{Proposal mechanism}
\begin{algorithm}[H]
\caption{Proposal mechanism}
\label{SA:Prop}
\begin{algorithmic}
\Procedure{ProposalMechanism}{} 
  \State $i \gets $ \emph{UniformRandom}($1,|N|$)
  \While {$i = j$}
    \State $j \gets $ \emph{UniformRandom}($1,|N|$)
  \EndWhile  
  \State \textbf{return} \emph{InsertTourComponent}($s,i,j$)   
\EndProcedure
\end{algorithmic}
\end{algorithm}

The probabilistic proposal mechanism consists of a random uniform choice of a candidate solution in the 1-shift insert neighborhood of the current solution.

\subsubsection{Acceptance criterion}
\begin{algorithm}[H]
\caption{Acceptance Criterion}
\label{saTSPTW}
\begin{algorithmic}
\Procedure{Acceptance criterion}{$s,s',T$}
\If {$f'(s) < f'(s')$}
  \State \textbf{return} $1$
\Else
  \State \textbf{return} $e^{\frac{f'(s)-f'(s')}{T}}$
\EndIf
\EndProcedure
\end{algorithmic}
\end{algorithm}

The acceptance criterion consists of an implementation of the metropolis condition, using a modified version of the evaluation function in order to take into account the number of constraint violations, while deciding whether to accept or not a worsening solution.

\begin{equation} \label{eq:metropolis}
  P(s,s',T) = \begin{cases}
               1 & f'(s) \le f'(s') \\
               e^{\frac{f'(s)-f'(s')}{T}} & \text{otherwise}
              \end{cases}
\end{equation}

where:

\begin{equation} \label{eq:modev}
  f'(s) = f(s) + \lambda \cdot \Omega(s)
\end{equation}

As \ref{eq:metropolis} show, downhill transitions are accepted with probability $1$, while the probability of accepting uphill transitions depends on the current temperature value $T$ and the difference between the values of modified evaluation function (\ref{eq:modev}) for the current ($s$) and proposed solution ($s'$).


\subsubsection{Temperature update}
\begin{algorithm}[H]
\caption{Update temperature according to annealing schedule}
\label{saTSPTW}
\begin{algorithmic}
\Procedure{Update}{$T,l$}
  \If{iteration \% $l$}
    \State \textbf{return} $\alpha \cdot T$
  \Else
    \State \textbf{return} $T$
  \EndIf
\EndProcedure
\end{algorithmic}
\end{algorithm}

The annealing schedule reduces the temperature of a coefficient $\alpha$ every $l$ time steps. 
This schedule allows to favor exploration (i.e. higher likelihood of accepting worsening step) at the beginning of the execution while focusing on exploitation as the number of iterations increases.

\subsubsection{Implementation parameters}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
\textbf{Parameter} & \textbf{Default} & \textbf{Selected} \\ \hline 
$\alpha$ & 0.95 & 0.8 \\\hline 
$T_0$ & 100 & 50 \\ \hline 
$l$ & 30000 & 1000 \\ \hline 
$LUBT$ & 100 & 100 \\ \hline 
$t_{\max}$ & 10[s] & 100[s] \\
\hline
\end{tabular}
\captionof{table}{SA - Algorithm parameters overview}
\label{saParameters}
\end{center}

The first parameters setting was determined according to the parameters described in \cite{ohlmann2007compressed}, also by computing the initial temperature as a function of the initial acceptance ratio.

Then, the proposed values has been determined by manually tuning the parameter to improved the quality of the solutions found by the algorithm.

I found that a cooling schedule with frequent temperature updates (smaller $l$), greater temperature reduction at each update (smaller $\alpha$) and lower initial temperature (smaller $T_0$), generally allows to obtain solution having a good quality even though the reduced temperature will sometimes cause the exploitation of sub-optimal solutions.

The maximum run-time $t_{\max}$ has been chosen equal to 100 seconds in order to analyze the run-time distribution of the algorithm over a sufficiently long run.
 
\end{homeworkProblem}		
